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Abstract. Building a training system for medical personnel to properly interpret fetal heart rate tracing requires 
developing accurate models that can relate various signal patterns to certain pathologies. In addition to modeling 
the fetal heart rate signal itself, the change of uterine pressure that bears strong relation to fetal heart rate and 
provides indications of maternal and fetal status should also be considered. In this work, we have developed a 
group of parametric models to simulate uterine contractions during labor and delivery. Through analysis of real 
patient records, we propose to model uterine contraction signals by three major components: regular contractions, 
impulsive noise caused by fetal movements, and low amplitude noise invoked by maternal breathing and measuring 
apparatus. The regular contractions are modeled by an asymmetric generalized Gaussian function and least 
squares estimation is used to compute the parameter values of the asymmetric generalized Gaussian function 
based on uterine contractions of real patients. Regular contractions are detected based on thresholding and 
derivative analysis of uterine contractions. Impulsive noise caused by fetal movements and low amplitude noise by 
maternal breathing and measuring apparatus are modeled by rational polynomial functions and Perlin noise, 
respectively. Experiment results show the synthesized uterine contractions can mimic the real uterine contractions 
realistically, demonstrating the effectiveness of the proposed algorithm. 


1. Introduction 

Uterine contractions are stimulated by uterine muscle 
cells. Uterine contraction (UC) variations reflect the 
physiological changes of the uterus during both 
pregnancy and labor [1], As a critical component in 
fetal heart rate (FHR) monitoring during labor and 
delivery, uterine contractions provide important 
information regarding maternal and fetal wellbeing. 
There are three main methods to record uterine 
contractions [2]: tocography, electrohysterography, 
and using intrauterine pressure catheter. 
Tocography measures the strain exerted by uterus on 
the maternal abdomen via external a tocotransducer. 
Electrohysterography records the electrical uterine 
activities from the maternal abdomen. The 
intrauterine pressure catheter (IUPC) measures 
intrauterine pressure invasively, and is mostly used 
during labor. Regardless of their differences, all 
three methods aim at providing records of the 
contraction patterns and their relationship to FHR. 
Proper interpretations of fetal heart rate and uterine 
contractions require special training, while monitoring 
of both are only available when pregnant women are 
hospitalized in parturiency. Building a training 
system that can simulate fetal heart rate and uterine 
contractions can help medical personnel learn critical 
patterns of both signals without putting patient in 
danger. To gain better understanding of uterine 
activities, it is necessary to develop mathematical 
models to quantitatively describe various uterine 
contraction patterns and this is the problem to be 
addressed in this paper. 

Even though FHR monitoring is now the standard 


practice during labor and delivery, surprisingly, there 
were only a few studies on uterine contraction 
modeling and simulation. Young used polynomials 
to model five characteristics of uterine contractions: 1) 
gradual onset, 2) a linear rising segment, 3) a plateau 
region, 4) a symmetrical fall, and 5) gradual offset, 
and fitted the simulated contractions with recorded 
IUPC data [3]. Their results matched their physical 
analysis. However, their simulations were not 
perfect especially in the tail region of the contraction 
curve. Vauge et al. [4] developed a system of 
differential equations that describe the dynamics of 
uterine pressure during human parturition. This 
method was based on three simplified assumptions: 
identical contractile properties of all myometrial cells, 
intrauterine pressure proportional to the number of 
contracted myometrial cells, and that all cells have 
three states, namely, contraction, recovery, and 
resting. Their model was simple and effective 

especially for normal contractions which begin in the 
fundus, reach the apex, and then proceed 
symmetrically downward toward the fundus. But this 
method did not consider the fact that asymmetry can 
occur when the uterine cells function independently 
causing ineffective uterine contractions and minimal 
dilatation [2]. Recently, Kemal et al. [5] employed 
two methods to simulate uterine contractions. The 
first one is based on the same mathematical model 
proposed by Vauge [4], The second approach is 
based on recorded patient data. They first applied 
Hilbert-Huang Transform (HHT) [6] to identify the 
contraction locations from real patient data, and then 
developed spatial-temporal simulations of uterine 
contractions. All these methods discussed above 
are capable of illustrating the dynamics of uterine 
activities. However, these methodologies are 
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deterministic and no comprehensive parameter 
estimation for these models was developed. 
Furthermore, the noise caused by fetal movements 
and maternal breathing were not considered and 
modeled. To address the problems in existing 
uterine modeling and simulation, the paper proposes 
a novel algorithm that integrates three major 
components: asymmetric generalized Gaussian 

function (AGGF) for modeling contractions, Perlin 
noise for modeling maternal breathing and instrument 
noise, and impulsive noise for modeling fetal 
movements. The parameters of the asymmetric 
generalized Gaussian function are estimated using 
the least square method based on detected uterine 
contractions from real patient records. 

The remainder of this paper is structured as follows. 
Section 2 first introduces the proposed asymmetric 
generalized Gaussian functions for modeling uterine 
contractions and then estimate the parameters of 
generalized Gaussian functions based detected 
contractions. Section 3 describes Perlin noise 
generation for maternal breathing and low amplitude 
noise modeling. Section 4 presents impulsive noise 
generation for simulating disturbances caused by 
fetal movements. Section 5 summarizes the 
simulation procedure and compares the simulation 
results. Section 6 concludes this paper and 
discusses future research directions. 

2. Uterine Contraction Modeling 

2.1 Asymmetric Generalized Gaussian Function 

A typical uterine contraction curve of a real patient is 
depicted as a continuous waveform in Figure 1 . This 
curve is characterized by a basal tone varying from 0 
to 20 units, and a deflection of the contraction curve 
above the baseline, whose amplitude and duration 
are within a certain range of values [2]. Moreover, 
one should note this curve is asymmetric. 



Figure 1: Typical uterine contractions. Also shown 
are parameters used for contraction detection. 

Since the asymmetry of contraction curve matches 


real cases of ineffective contractions and comprises 
the symmetrical case, we propose to use an 
asymmetric generalized Gaussian function to model 
uterine contractions as follows. 

fit) = A, exp |[« {t-t,)-u (t - 1 0 )] 
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where the parameters of (1) are shown below. 

A/, A r Amplitudes for the left and right sides 

a/, a r Exponents for the left and right sides 

Pi, p r Variances for the left and right sides 
u{t) A unit step function 
h t r Left and right cut off time 
t 0 The position where f(t) reaches its maximum 

b(t ) Baseline representing some basal strain 

exerted by the uterine muscle when 
contractions do not occur 

To simulate uterine contractions, we need to 
determine the range within which the above 
parameters lie and how they vary with time. So in 
the next step we will detect uterine contractions from 
real patient data and estimate the parameters for the 
asymmetric generalized Gaussian function from the 
detection results. 

2.2 Uterine Contraction Detection 

Several methods were proposed to detect uterine 
contractions for different purposes. Radhakrishnan 
et al. [7] developed a higher-order zero crossing 
based method and studied the frequency of 
occurrence of contractions in different pregnancy 
stages. Novak et al. [8] described two UC detection 
approaches: amplitude- and derivative-based 

algorithm. By comparing the results from these two 
methods, they suggested combining both methods 
together to achieve better detection results. Aiming 
at quantitatively analyzing uterine contractions in time 
domain, Jezewski et al. [9] introduced a statistical 
method to determine the threshold and also 
considered duration condition for UC detection. 
These methods were used to calculate the regular 
parameters of uterine contractions such as amplitude, 
duration, frequency of occurrence. Our method is 
based on a combination of the last two methods, in 
which thresholding is first performed to detect the 
presence of uterine contractions, and 
derivative-based method is applied subsequently to 
include the samples whose amplitudes are below the 
threshold, but still belonging to the contractions. 

The algorithm proposed by Jezewski et al. [9] starts 
with low pass filtering with cutoff frequency of 0.04 Hz 
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to suppress the artifacts caused by fetal movements 
and maternal breathing. Then the record is 

analyzed by using a moving window with a length of 4 
minutes and 1 -minute step. Within each window, the 
histogram of uterine pressure samples is constructed 
first and the mode of the histogram is then selected 
as the baseline value. Finally, the threshold level is 
set as 10 units above the baseline and the validity of 
data segment is examined. A valid contraction 
should remains above the threshold level for a 
duration longer than 30 seconds and the amplitude of 
contraction exceeds 20 units. 


(a) Rj» UC (UU 



Figure 2: Intermediate results in uterine contraction 
detection, (a) Raw UC data, (b) Filtered UC. (c) 
UC detection by thresholding, (d) The derivatives of 
filtered UC. (e) Improved UC detection results, (f) 
The residual between raw UC and filtered UC which is 
detected. 

After the uterine contractions are detected, there are 
still two remaining problems. The first is that the 
amplitude of detected samples are mostly above the 
threshold, thus the estimated parameters may not be 
accurate enough to model the tail areas. We need to 
include more samples from tail areas. On the other 
hand, the tail areas are prone to other various 
sources of noise. Before including them for parameter 
estimation, we must differentiate contaminated 
samples and uncontaminated or less-contaminated 
samples. Since the derivatives in noisy area change 
drastically, we propose to employ derivative-based 
detection method to differentiate them. In other 
words, for less contaminated contractions, the 
derivative of the left side of the contraction should be 
positive, while the derivative of right side should be 


negative. Starting from the peak of the contraction, 
we proceed to its left and compute the derivatives of 
the smoothed uterine pressure signal and then search 
the position where the derivative changed to negative 
for left side curve, and denote it as the tail point. The 
sample points between the peak and the tail point will 
be utilized for curve parameters estimation. The 
same principle applies to the right side of contraction 
with positive derives being searched. 

The results of uterine contraction detection are shown 
Figure 2, where (a) is the original uterine pressure 
signal; (b) is its filtered version, the red solid line 
represents the base line; (c) illustrates the sample 
whose values are above the threshold; (d) is the 
derivatives of filtered UC data; (e) is the improved 
uterine contraction detection result; and (f) shows the 
residual between the detected contractions and 
original signal (a). 


2.3 Uterine Contraction Parameter Estimation 

After the uterine contractions are detected, we need 
to compute the parameters of the asymmetric 
generalized Gaussian function (1) in order to simulate 
the detected contractions. Here we only consider 
the left half of the asymmetric generalized Gaussian 
function (1), that is, 

/(/) = 4exp|-l^-^L| + 6 < [«(/-/,)-«(/-/ 0 )] 

( 2 ) 

for the parameter estimation, while the right half can 
be handled similarly and thus is omitted in this paper. 


First, the baseline b t is estimated as the contraction 
value at the onset point t t which is the lowest point 
between two filtered contractions. Thus 


4 =/(' 0 )-V 


After simple manipulation, the remaining known 
variables and unknown parameters are separated by 
taking the logarithm of both sides of (2), 
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(3) 


Since < A, , adding minus sign to both sides 

of (3) and taking logarithm again gives 
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Now we denote the detected uterine contraction data 
set as j, \< i < N , where N is the 

number of samples. Substituting the data set into 
equation (4) and writing each term in matrix form, we 
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Then equation (4) can be transformed into a concise 
linear equation represented by a matrix product, 

Y = Xp, (6) 

where X and Y are matrices containing detected 
sample information, p is the parameter vector to be 
estimated. The parameter estimation problem can 
be expressed as the following minimization problem, 

P„ =argmin||Y-Xp||J. (7) 

peR 2 


where p o is the solution. This is a least square 
estimation problem whose solution is 


P 0 = (X ? X) ' X 7 Y . (8) 


Thus we have obtained the parameter estimator for 
the proposed asymmetric generalized Gaussian 
function for modeling uterine contractions. 


3. Perlin Noise 

The low amplitude noise caused by maternal 
breathing and measuring apparatus are random yet 
exhibiting both low and high frequency characteristics 
as shown in Figures 2(a) and (f). Common random 
number generators cannot be used directly to 
generate such noise, since they are too random to 
exhibit the natural outlook of continuity and 
self-similarity of the noise. To address this problem, 
we propose to use Perlin noise generator [10] to 
simulate the low amplitude noise. 

A Perlin noise generator is composed of two 
components: a noise function and an interpolation 
function. The basic procedure of Perlin noise 
generation is 

(1) Generate a series of random numbers A? y of 
length N from uniform distribution U[-0.1. 0.1]; 

(2) Decimate the series to size N/2\ 

(3) Upsample the decimated series to size N by 
B-spline interpolation and increase the amplitude 
of the new series by a factor of p; 

(4) Repeat (2) and (3) until reaching the specified 

level s, and we obtain a set of series { n L n 2 ,...,n s }. 


By adding the above generated series together, we 
have the Perlin noise 

4(0=Zp'«,(0* ( 9 ) 

i=0 

where fe[l, N], and p is called persistence factor 
which controls the amplitude. 





(c) 



Figure 3: (a) - (f) are waveforms corresponding to 6 
levels noise series {n, n 2 .....n 6 ). (g) is the Perlin noise. 
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One fact worth of explaining is the mechanism for 
how the continuous noise is produced. Since n, is 
generated from i.i.d. (independent and 
identically-distributed) uniform distribution, it is 
actually white noise. Through repeated 
downsampling and upsampling, the new noise series 
become the low-pass filtered version of the previous 
noise series. In other words, decimation by a factor 
of 2 reduces one half the Nyquist frequency of 
previous noise, then B-spline interpolation restores 
the sample number of noise without incurring new 
frequency contents. Thus the noise n s generated in 
the last step, occupies the lowest frequency band. 
So the waveforms of the noise series from n x to n s , 
become increasingly smoother. Figure 3 illustrates 
the components of the Perlin noise generated in this 
work, in which (a) is the noise generated from uniform 
distribution U[-0.1, 0.1], (b) - (f) are low-pass filtered 
version of noise generated from one level before with 
p = 2. Note that the sub-figures in Figure 3 have 
different vertical scales and low-frequency 
components have much larger amplitudes than 
high-frequency components. The final synthesized 
Perlin noise is shown in Figure 3(g). 

4. Impulsive Noise Modeling 

The last component to simulate is the spikes or large 
magnitude impulsive noise in uterine contractions, 
which suggest possible fetal movements [2]. After 
examining the shape of those spikes, we propose to 
use the following rational polynomial function to 
model them 


m= 


ft-t v 
1 *0 


1 + 


V o j 


( 8 ) 


where a is the amplitude which follows uniform 
distribution with value between 0 and 50 units, b is the 
scale parameter following uniform distribution with 
value between 0 and 10, and t 0 is the spike position. 
A typical simulation result of the impulsive noise is 
shown in Figure 4. 



Figure 4: Impulsive noise simulation 


5. Uterine Contraction Synthesis 

Uterine contraction simulation is finalized by 
superimposing the components generated by the 
asymmetric generalized Gaussian function, Perlin 
noise, and impulsive noise. First, a typical segment 
of 20 minutes was extracted from real patient record 
and the parameters for the asymmetric generalized 
Gaussian function were estimated. Figure 5 
compares the parameters estimated from two 
detected data sets. The first data set contains only 
samples whose amplitude is above the threshold. 
The second data set is the expanded version of the 
first set by incorporating samples whose amplitude 
are below the threshold but are belong to the 
contraction. Figures 5 (a) and (b) plot parameters 
estimated from the first set. Figures 5(c) and (d) are 
estimated from the second set. It can be seen that 
the parameters of the left side and right side are 
different, validating the asymmetry of uterine 
contractions. 

To compare the impact of different parameters on 
uterine contraction simulation, we select one set of 
parameters from Figures 5(a) and (b), and another 
set from Figures 5(c) and (d) to simulate two 
contractions, namely, contraction 1 (green dotted line) 
and contraction 2 (red dashed line), as shown in 
Figure 5(e). It can be seen that the red dash line 
achieves a better fit to the real contraction at both the 
peak and tail area of the contraction and has smaller 
normalized root mean square error (NRMSE), which 
is defined as 


NRMSE = 


j Zlto »-fM 
US <AO)>-gta(/.(*)}' 


( 10 ) 


This result indicates that the derivative-based 
analysis could help to recruit more effective data for 
parameter estimation, and thus increase the 
estimation accuracy. Finally, by adding together the 
asymmetric Gaussian function, the Perlin noise, and 
the impulsive noise, we obtain the final simulation 
results in Figure 6(b). Comparing the original uterine 
contraction record in Figure 6(a) and the simulated 
results in Figure 6(b), it can be easily seen that the 
proposed algorithm is very effective and produces 
superb results. 
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Figure 5: Parameters estimation and uterine 

contraction simulation, (a) is the scatter plot of a, 
and (b) is the scatter plot of p for left and right curves, 
(c) and (d) are a and p plots of the second data set. 
(e) shows the simulations of one uterine contraction 
from two data sets. 


(a) recorded data 



(b) simulation results 



Figure 6: Simulation results, (a) The original uterine 
contraction recording, (b) Simulated uterine 
contraction based on (a). It can be seen that the 
proposed algorithm is very effective, producing 
realistic simulations. 

6. Conclusion 

This paper proposed a set of parametric models to 
simulate uterine contractions. The proposed 
algorithm contains three major components: AGGF 
model for contractions, Perlin noise for maternal 
breathing and instrument noise, and rational 
polynomial functions for fetal movements. 


Parameter estimation for the asymmetric generalized 
Gaussian function is derived based on real uterine 
contractions. The proposed algorithm is effective 
and produces realistic simulation results. Future 
work includes parameter estimation for the Perlin 
noise and impulsive noise in order to further improve 
the quality of uterine contraction simulations. 
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